%%  Assignment 2U8 for WI3097TU-p
% This is the published matlab-file

% Authors:      Martin Tichem, Bas des Tombe
% Date:         2012-06-04 
% Available on: https://code.google.com/p/wi3097tu/

% Clear Memory for Initialization
clear all;close all;clc

%% Question 9

clear all;
fprintf('\n___________________\n### QUESTION 9 ###\n')

% Define constants
const.k_d   = 0.24;
const.k_s   = 0.19;
const.mu    = 4;
const.Pk    = 8;
const.q     = 0.01;
const.P_in  = .5;

h       = 4;                    % Stapgrootte
tmax    = 40;
t       = 0:h:tmax;             % Tijdsdomein
winit   = [3/80;1/16;9/200];    % Initial values

% Function call
w=funcEM(const,h,tmax,winit);

% Output
% option = input('What would you like as an output?\n Option 1: Fancy UI-Table\n Option 2: Horizontal Table\n Option 3: Vertical Table\n');
option = 3;

switch option
    case 1
        % Mooie tabel
        f = figure('Position',[200 200 1040 150]);
        rnames = {'Tijd [uren]','Algen [mol/m3]','Fosfaten [mol/m3]','Detritus [mol/m3]'};
        table = uitable('Parent',f,'Data',single([t;w]),'RowName',rnames,'ColumnWidth','auto','Position',[20 20 1020 100]);

    case 2
        % Minder mooie tabel
        fprintf('Tijd [uren]')
        for n = 1:length(t)
            fprintf('\t\t\t%6.0f',t(n));
        end

        fprintf('\n Algen [mol/m3]\t')
        for n = 1:length(t)
            fprintf('\t%1.7e',w(1,n));
        end
        fprintf('\n Fos\n\t faten [mol/m3]')
        for n = 1:length(t)
            fprintf('\t%1.7e',w(2,n));
        end
        fprintf('\n Detritus [mol/m3]')
        for n = 1:length(t)
            fprintf('\t %1.7e',w(3,n));
        end
        fprintf('\n')% 
    
    case 3
        % Tabel
        fprintf('\n---TABLE---\n    x   Algen\t\tFosfaten\t Detritus\n');
        for n = 1:length(t)
            fprintf('%6.0f \t%1.7e \t%1.7e \t %1.7e \n',t(n),w(1,n),w(2,n),w(3,n));
        end
        
    otherwise
        error('Invalid input');
end
 
%%  Question 10
clear all;
fprintf('\n_______\n\t ____________\n### QUESTION 10 ###\n\n')

% Define constants
const.k_d   = 0.24;
const.k_s   = 0.19;
const.mu    = 4;
const.Pk    = 8;
const.q     = 0.01;
const.P_in  = .5;

h       = 1;                    % Stapgrootte
tmax    = 250;
t       = 0:h:tmax;             % Tijdsdomein
winit   = [3/80;1/16;9/200];    % Initial values

% Function call
w=funcEM(const,h,tmax,winit);

% Plot
g       = figure();
plot(t,w(1,:),t,w(2,:),t,w(3,:));
legend('Algen','Fosfaten','Detritus');
xlabel('Tijd [uren]')
ylabel('[mol/m^3]')

fprintf('See figure 1\n');
%%  Question 11

clear all;
fprintf('\n___________________\n### QUESTION 11 ###\n')

% Define constants
const.k_d   = 0.24;
const.k_s   = 0.19;
const.mu    = 4;
const.Pk    = 8;
const.q     = 0.01;
const.P_in  = .5;

h_1     = 1;                    % Stapgrootte 1
h_2     = 2;                    % Stapgrootte 2
tmax    = 50;
t_1     = 0:h_1:tmax;           % Tijdsdomein 1
t_2     = 0:h_2:tmax;           % Tijdsdomein 2
winit	= [3/80;1/16;9/200];	% Initial values

% Function call; Calculate error
wl1=funcEM(const,h_1,tmax,winit);
wl2=funcEM(const,h_2,tmax,winit);

p=2;
alpha=1/(1-2^p);
error=alpha*(wl2(:,end)-wl1(:,end))

%%  Question 12 

clear all;
fprintf('___________________\n### QUESTION 12 ###\n\n')

% Define constants
const.k_d   = 0.24;
const.k_s   = 0.19;
const.mu    = 4;
const.Pk    = 8;
const.q     = 0.01;

h       = 1;                    % Stapgrootte
tmax    = 1000;
t       = 0:h:tmax;             % Tijdsdomein
winit   = [3/80;1/16;9/200];    % Initial values

P2                          = .1;
const.P_in                  = zeros(1,length(t)); % Allocating P_in
const.P_in(1:100)           = 0.5;
const.P_in(101:400)         = linspace(0.5,P2,300/h);
const.P_in(401:length(t))   = P2;

[w,itilde]                  = funcEM2(const,h,tmax,winit);
fprintf('The value of t-tilde is: %i [h]\n',t(itilde));

%%  Question 13

clear all;
fprintf('\n___________________\n### QUESTION 13 ###\n\n')

% Define constants
const.k_d   = 0.24;
const.k_s   = 0.19;
const.mu    = 4;
const.Pk    = 8;
const.q     = 0.01;

h       = 1;                    % Stapgrootte
tmax    = 1000;
t       = 0:h:tmax;             % Tijdsdomein
winit   = [3/80;1/16;9/200];    % Initial values

P2      = [.05, .06, .1, .15, .2];

for i = 1:length(P2);
    const.P_in                  = zeros(1,length(t)); % Allocating P_in
    const.P_in(1:100)           = 0.5;
    const.P_in(101:400)         = linspace(0.5,P2(i),300/h);
    const.P_in(401:length(t))   = P2(i);
    
    [w,itilde]                  = funcEM2(const,h,tmax,winit);
    if itilde ~= 1;
        fprintf('For P2 is %1.3f, the value of t-tilde is %i hours\n',P2(i),t(itilde));
    else
        fprintf('For P2 is %1.3f, the algae concentration does not drop below .005 mol/m3\n\n',P2(i));
    end
end